#include "projectionlambert.h"
//数学库
#include<cmath>
projectionLambert::projectionLambert()
{

}

projectionLambert::projectionLambert(double b0,double b1,double l0,double l1)
{
    this->B0=b0;//最低纬度
    this->B1=b1;//最高纬度
    this->L0=l0;//最小经度
    this->L1=l1;//最大经度

    //各种计算用常数
    double alpha=1/298.257;
    a=6378137;
    e=sqrt(2*alpha-pow(alpha,2));

    Lm=(L0+L1)/2;
    stdB0=B0+0.16*(B1-B0);
    stdB1=B1-0.12*(B1-B0);

    r0=(a*cos(stdB0*PI/180))/(pow((1-(e*e*sin(stdB0*PI/180)*sin(stdB0*PI/180))),0.5));
    r1=(a*cos(stdB1*PI/180))/(pow((1-(e*e*sin(stdB1*PI/180)*sin(stdB1*PI/180))),0.5));
    rs=(a*cos(B0*PI/180))/(pow((1-(e*e*sin(B0*PI/180)*sin(B0*PI/180))),0.5));

    u0=tan((PI/4)+(stdB0/2*PI/180))*pow((1-e*sin(stdB0*PI/180))/(1+e*sin(stdB0*PI/180)),e/2);
    u1=tan((PI/4)+(stdB1/2*PI/180))*pow((1-e*sin(stdB1*PI/180))/(1+e*sin(stdB1*PI/180)),e/2);
    us=tan((PI/4)+(B0/2*PI/180))*pow((1-e*sin(B0*PI/180))/(1+e*sin(B0*PI/180)),e/2);

    alphaC=(log(r0)-log(r1))/(log(u1)-log(u0));
    C=r1*pow(u1,alphaC)/alphaC;

    rhoS=rs*pow(us,alphaC)-alphaC*C;
}

projectionLambert::~projectionLambert()
{

}

void projectionLambert::getXY(double L, double B, double &x, double &y)
{
    //L: 经度, B: 纬度
    double rho,delta,u,l;
    double mul=1;
    l=(L/mul-Lm)*PI/180;
    u=tan((PI/4)+(B/2*PI/mul/180))*pow((1-e*sin(B*PI/mul/180))/(1+e*sin(B*PI/mul/180)),e/2);
    rho=C/pow(u,alphaC);
    delta=alphaC*l;

    y=rhoS-rho*cos(delta);
    x=rho*sin(delta);
}
